This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
library(raster)
size_dat=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/fire_growth_5days_v4.txt", header=T, row.names=NULL)
size_dat=as.data.frame(size_dat)
size_dat=size_dat[-1,]
size_dat=size_dat[,-1]
colnames(size_dat)=c("firename","year","cause","size1","size2","size3","size4","size5","final_firesize","mean_precip1","mean_precip2","mean_precip3","mean_precip4","mean_precip5","mean_tmax1","mean_tmax2","mean_tmax3","mean_tmax4","mean_tmax5","mean_tmean1","mean_tmean2","mean_tmean3","mean_tmean4","mean_tmean5","mean_vpdmax1","mean_vpdmax2","mean_vpdmax3","mean_vpdmax4","mean_vpdmax5","mean_windspeed1","mean_windspeed2","mean_windspeed3","mean_windspeed4","mean_windspeed5","landcover","ecosystem","biomass","elevation")
size_dat2 <- data.frame(lapply(size_dat, function(x) as.numeric(as.character(x))))
NAs introduced by coercion
size_dat2$human = 0
size_dat2$human[size_dat2$cause !=1 & size_dat2$cause !=14 & size_dat2$cause !=17]=1
size_dat2$human[size_dat2$cause ==1 ]=2
pro1 =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1 ), ]
pro2 =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1 ), ]
t.test(pro1$size1,pro2$size1)
Welch Two Sample t-test
data: pro1$size1 and pro2$size1
t = -2.0397, df = 74.381, p-value = 0.04493
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-22.0986230 -0.2593373
sample estimates:
mean of x mean of y
3.056385 14.235365
t.test(pro1$size2,pro2$size2)
Welch Two Sample t-test
data: pro1$size2 and pro2$size2
t = -2.758, df = 73.138, p-value = 0.007341
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-32.702045 -5.266208
sample estimates:
mean of x mean of y
8.254394 27.238521
t.test(pro1$size3,pro2$size3)
Welch Two Sample t-test
data: pro1$size3 and pro2$size3
t = -2.9002, df = 58.203, p-value = 0.005254
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-51.965993 -9.526741
sample estimates:
mean of x mean of y
14.10484 44.85121
t.test(pro1$size4,pro2$size4)
Welch Two Sample t-test
data: pro1$size4 and pro2$size4
t = -3.1797, df = 53.078, p-value = 0.002461
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-74.99798 -16.98040
sample estimates:
mean of x mean of y
17.68430 63.67349
t.test(pro1$size5,pro2$size5)
Welch Two Sample t-test
data: pro1$size5 and pro2$size5
t = -3.4218, df = 39.443, p-value = 0.001462
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-117.71273 -30.26844
sample estimates:
mean of x mean of y
19.71795 93.70853
pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1), ]
length(pro$year)
[1] 68
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1), ]
length(pro$year)
[1] 77
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 2), ]
length(pro$year)
[1] 18
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 2), ]
length(pro$year)
[1] 9
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
length(!is.na(pro$size5))
[1] 9
pro =size_dat2[which(size_dat2$human == 1 & size_dat2$ecosystem == 6), ]
length(pro$year)
[1] 41
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
pro =size_dat2[which(size_dat2$human == 2& size_dat2$ecosystem == 6), ]
length(pro$year)
[1] 79
boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)
NA
NA
plot fire size map QGIS
library(data.table)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
data.table 1.13.0 using 1 threads (see ?getDTthreads). Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********
Attaching package: ‘data.table’
The following object is masked from ‘package:raster’:
shift
DT= data.table(res)
fire_size = DT[ , .SD[which.min(growth_km)], by = firename]
Error in .checkTypos(e, names_x) :
Object 'growth_km' not found amongst bi, erc, etr, fm100, fm1000 and 23 more
library(raster)
dr1 =shapefile("/Users/stijnhantson/Documents/data/FRAP/fire18_1.shp")
dr1$YEAR_=as.numeric(as.character(dr1$YEAR_))
dr1$Shape_Area=as.numeric(as.character(dr1$Shape_Area))
dr1$ALARM_DATE = as.Date(dr1$ALARM_DATE )
dr1=dr1[!is.na(dr1$YEAR_), ]
dr1=dr1[dr1$YEAR_>2011,]
dr2=dr1[dr1$YEAR_ == 2013,]
plot(table(dr2$ALARM_DATE), xlim=c(as.Date("2013/01/01"),as.Date("2013/12/31")))
table(factor(dr2$ALARM_DATE, levels = as.Date("2013/01/01"):as.Date("2013/12/31")))
15706 15707 15708 15709 15710 15711 15712 15713 15714 15715 15716 15717 15718 15719 15720
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15721 15722 15723 15724 15725 15726 15727 15728 15729 15730 15731 15732 15733 15734 15735
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15736 15737 15738 15739 15740 15741 15742 15743 15744 15745 15746 15747 15748 15749 15750
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15751 15752 15753 15754 15755 15756 15757 15758 15759 15760 15761 15762 15763 15764 15765
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15776 15777 15778 15779 15780
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15781 15782 15783 15784 15785 15786 15787 15788 15789 15790 15791 15792 15793 15794 15795
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15806 15807 15808 15809 15810
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15811 15812 15813 15814 15815 15816 15817 15818 15819 15820 15821 15822 15823 15824 15825
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15826 15827 15828 15829 15830 15831 15832 15833 15834 15835 15836 15837 15838 15839 15840
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15841 15842 15843 15844 15845 15846 15847 15848 15849 15850 15851 15852 15853 15854 15855
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15856 15857 15858 15859 15860 15861 15862 15863 15864 15865 15866 15867 15868 15869 15870
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15871 15872 15873 15874 15875 15876 15877 15878 15879 15880 15881 15882 15883 15884 15885
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15886 15887 15888 15889 15890 15891 15892 15893 15894 15895 15896 15897 15898 15899 15900
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15901 15902 15903 15904 15905 15906 15907 15908 15909 15910 15911 15912 15913 15914 15915
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15916 15917 15918 15919 15920 15921 15922 15923 15924 15925 15926 15927 15928 15929 15930
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15931 15932 15933 15934 15935 15936 15937 15938 15939 15940 15941 15942 15943 15944 15945
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15946 15947 15948 15949 15950 15951 15952 15953 15954 15955 15956 15957 15958 15959 15960
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15961 15962 15963 15964 15965 15966 15967 15968 15969 15970 15971 15972 15973 15974 15975
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15976 15977 15978 15979 15980 15981 15982 15983 15984 15985 15986 15987 15988 15989 15990
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15991 15992 15993 15994 15995 15996 15997 15998 15999 16000 16001 16002 16003 16004 16005
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
16006 16007 16008 16009 16010 16011 16012 16013 16014 16015 16016 16017 16018 16019 16020
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
16021 16022 16023 16024 16025 16026 16027 16028 16029 16030 16031 16032 16033 16034 16035
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
16036 16037 16038 16039 16040 16041 16042 16043 16044 16045 16046 16047 16048 16049 16050
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
16051 16052 16053 16054 16055 16056 16057 16058 16059 16060 16061 16062 16063 16064 16065
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
16066 16067 16068 16069 16070
0 0 0 0 0
human lightning unknown ratio through time
library(raster)
library(rgdal)
package ‘rgdal’ was built under R version 3.6.2rgdal: version: 1.5-12, (SVN revision 1018)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE
Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
Linking to sp version:1.4-2
Overwritten PROJ_LIB was /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
library(rgeos)
package ‘rgeos’ was built under R version 3.6.2rgeos version: 0.5-3, (SVN revision 634)
GEOS runtime version: 3.7.2-CAPI-1.11.2
Linking to sp version: 1.4-1
Polygon checking: TRUE
library(sf)
frap=readOGR("/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp")
OGR data source with driver: ESRI Shapefile
Source: "/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp", layer: "fire19_1"
with 20820 features
It has 18 fields
Integer64 fields read as strings: CAUSE C_METHOD OBJECTIVE
shape = shapefile("/Users/stijnhantson/Documents/data/veg_california/ca_eco_l3/ca_eco_l3.shp")
shape = spTransform(shape,crs(frap))
frap = st_make_valid(st_read("/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp"))
Reading layer `fire19_1' from data source `/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp' using driver `ESRI Shapefile'
GDAL Message 1: organizePolygons() received an unexpected geometry. Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry. Return arguments as a collection.GDAL Message 1: Geometry of polygon of fid 19691 cannot be translated to Simple Geometry. All polygons will be contained in a multipolygon.
Simple feature collection with 20820 features and 18 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -373237.5 ymin: -604727.6 xmax: 519987.8 ymax: 518283.7
CRS: 3310
types <- vapply(sf::st_geometry(frap), function(x) {
class(x)[2]
}, "")
polys <- frap[ grepl("*POLYGON", types), ]
frap=sf:::as_Spatial(polys)
frap=gBuffer(frap, byid=T, width=0.0)
eco_inter=intersect(frap,shape)
NA
library(raster)
#library(rgdal)
daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/final_dataset_V5.txt",header=T)
res=as.data.frame(daily_res)
res$mean_ros =as.numeric(as.character(res$mean_ros))
res$max_ros =as.numeric(as.character(res$max_ros))
res$median95_ros =as.numeric(as.character(res$median95_ros))
res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
#res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$growth =as.numeric(as.character(res$growth))
res$total_area =as.numeric(as.character(res$total_area))
res$mean_frp =as.numeric(as.character(res$mean_frp))
res$frp_95 =as.numeric(as.character(res$frp_95))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))
res$biomass =as.numeric(as.character(res$biomass))
res$year =as.numeric(as.character(res$year))
res$month =as.numeric(as.character(res$month))
res$doy_out =as.numeric(as.character(res$doy_out))
res$mean_dnbr =as.numeric(as.character(res$mean_dnbr))
res$mean_rdnbr =as.numeric(as.character(res$mean_rdnbr))
res$mean_bas =as.numeric(as.character(res$mean_bas))
res$median_dnbr =as.numeric(as.character(res$median_dnbr))
res$median_rdnbr =as.numeric(as.character(res$median_rdnbr))
res$median_bas =as.numeric(as.character(res$median_bas))
res$q95_dnbr =as.numeric(as.character(res$q95_dnbr))
res$q95_rdnbr =as.numeric(as.character(res$q95_rdnbr))
res$q95_bas =as.numeric(as.character(res$q95_bas))
res = res[-1,]
res$per_ba = res$growth/res$total_area
res$growth_km =res$growth/1000000
res$human = 0
res$human[res$cause !=1 & res$cause !=14 & res$cause !=17]=1
res$human[res$cause ==1 ]=2
res$ros_km = (res$median95_ros*24)/1000
res$ros_mean_km = (res$mean_ros*24)/1000
is there relation between days untill 75% and fire size
dr1 =shapefile("/Users/stijnhantson/Documents/data/FRAP/fire18_1.shp")
dr1$YEAR_=as.numeric(as.character(dr1$YEAR_))
dr1$Shape_Area=as.numeric(as.character(dr1$Shape_Area))
dr1=dr1[!is.na(dr1$YEAR_), ]
dr1=dr1[dr1$YEAR_>2011,]
fi = dr1[dr1$OBJECTID==17495,]
dr1=subset(dr1, GIS_ACRES>300)
summary(dr1)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -332881.6 349023.3
y -601475.0 462306.3
Is projected: TRUE
proj4string :
[+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs]
Data attributes:
OBJECTID YEAR_ STATE AGENCY UNIT_ID FIRE_NAME INC_NUM
Min. :11970 Min. :2012 Length:515 Length:515 Length:515 Length:515 Length:515
1st Qu.:17777 1st Qu.:2014 Class :character Class :character Class :character Class :character Class :character
Median :19129 Median :2016 Mode :character Mode :character Mode :character Mode :character Mode :character
Mean :19125 Mean :2015
3rd Qu.:20512 3rd Qu.:2017
Max. :21127 Max. :2018
ALARM_DATE CONT_DATE CAUSE COMMENTS REPORT_AC GIS_ACRES
Length:515 Length:515 Length:515 Length:515 Min. : 0.0 Min. : 303.8
Class :character Class :character Class :character Class :character 1st Qu.: 504.1 1st Qu.: 612.4
Mode :character Mode :character Mode :character Mode :character Median : 1651.5 Median : 1707.6
Mean : 12021.3 Mean : 12138.1
3rd Qu.: 6970.0 3rd Qu.: 7068.0
Max. :410203.0 Max. :410202.5
NA's :67
C_METHOD OBJECTIVE FIRE_NUM Shape_Leng Shape_Area
Length:515 Length:515 Length:515 Min. : 4635 Min. :1.229e+06
Class :character Class :character Class :character 1st Qu.: 10017 1st Qu.:2.478e+06
Mode :character Mode :character Mode :character Median : 17475 Median :6.910e+06
Mean : 41493 Mean :4.912e+07
3rd Qu.: 37905 3rd Qu.:2.860e+07
Max. :445282 Max. :1.660e+09
res75 = res[res$per_ba > 0.999,]
peak_day1 = as.data.frame(aggregate(res75$fire_day, by = list(res75$firename,res75$year), min))
p=13
year1=0
first=0
last=0
len_d=0
day75=0
firenam1=0
name_fire=0
size=0
k=0
for (p in 1:(length(peak_day1$Group.1)))
{
print(p)
year = peak_day1[p,2]
firenam = as.character(peak_day1[p,1])
fir = dr1[which(dr1$YEAR_==year & dr1$FIRE_NAME==firenam),]
if (length(fir)>1){
maxsize = max(fir$GIS_ACRES)
fir = fir[fir$GIS_ACRES == maxsize,]
}
if (!is.na(fir$ALARM_DATE) | !is.na(fir$ALARM_DATE)){
k=k+1
day75[k] = peak_day1[p,3]
firenam1[k]=firenam
year1[k] = year
first[k]= fir$ALARM_DATE
last[k] = fir$CONT_DATE
len_d[k] = as.numeric(as.Date(fir$CONT_DATE) -as.Date(fir$ALARM_DATE))
name_fire[k]= fir$OBJECTID
size[k] = fir$GIS_ACRES
}
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
[1] 96
[1] 97
[1] 98
[1] 99
[1] 100
[1] 101
[1] 102
[1] 103
[1] 104
[1] 105
[1] 106
[1] 107
[1] 108
[1] 109
[1] 110
[1] 111
[1] 112
[1] 113
[1] 114
[1] 115
[1] 116
[1] 117
[1] 118
[1] 119
[1] 120
[1] 121
[1] 122
[1] 123
[1] 124
[1] 125
[1] 126
[1] 127
[1] 128
[1] 129
[1] 130
[1] 131
[1] 132
[1] 133
[1] 134
[1] 135
[1] 136
[1] 137
[1] 138
[1] 139
[1] 140
[1] 141
[1] 142
[1] 143
[1] 144
[1] 145
[1] 146
[1] 147
[1] 148
[1] 149
[1] 150
[1] 151
[1] 152
[1] 153
[1] 154
[1] 155
[1] 156
[1] 157
[1] 158
[1] 159
[1] 160
[1] 161
[1] 162
[1] 163
[1] 164
[1] 165
[1] 166
[1] 167
[1] 168
[1] 169
[1] 170
[1] 171
[1] 172
[1] 173
[1] 174
[1] 175
[1] 176
[1] 177
[1] 178
[1] 179
[1] 180
[1] 181
[1] 182
[1] 183
[1] 184
[1] 185
[1] 186
[1] 187
[1] 188
[1] 189
[1] 190
[1] 191
[1] 192
[1] 193
[1] 194
[1] 195
[1] 196
[1] 197
[1] 198
[1] 199
[1] 200
[1] 201
[1] 202
[1] 203
[1] 204
[1] 205
[1] 206
[1] 207
[1] 208
[1] 209
[1] 210
[1] 211
[1] 212
plot(len_d,day75, xlab="fire duration FRAP (days)", ylab="fire duration VIIRS (days)")
plot mean vs max fire rate-of-spread
#summary(res)
plot(res$ros_km,res$ros_mean_km, xlab="maximum fire-spread-rate (km/day",ylab="mean fire-spread-rate (km/day)")
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/mean_vs_max_ros_v1.pdf", width = 7, height = 7)
plot(res$ros_km,res$ros_mean_km, xlim=c(0,25),ylim=c(0,10), xlab="fire rate-of-spread (km/day)",ylab="mean fire rate-of-spread (km/day)", cex.lab=1.3,cex.axis = 1.25)
dev.off()
quartz_off_screen
2
difference between human and lightnign fires
me=0
me1=0
days = c("day1","day2","day3","day4","day5")
days1 = c("1","2","3","4","5")
pro1 = res[res$fire_day == 1 & res$human == 1,]
pro2 = res[res$fire_day == 2 & res$human == 1,]
pro3 = res[res$fire_day == 3 & res$human == 1,]
pro4 = res[res$fire_day == 4 & res$human == 1,]
pro5 = res[res$fire_day == 5 & res$human == 1,]
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
pro1h = res[res$fire_day == 1 & res$human == 2,]
pro2h = res[res$fire_day == 2 & res$human == 2,]
pro3h = res[res$fire_day == 3 & res$human == 2,]
pro4h = res[res$fire_day == 4 & res$human == 2,]
pro5h = res[res$fire_day == 5 & res$human == 2,]
me1[1] =mean(pro1h$growth_km,na.omit=T)
me1[2] =mean(pro2h$growth_km,na.omit=T)
me1[3] =mean(pro3h$growth_km,na.omit=T)
me1[4] =mean(pro4h$growth_km,na.omit=T)
me1[5] =mean(pro5h$growth_km,na.omit=T)
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure1_v2.pdf", width = 10, height = 5)
par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^2*')'),ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)
text(0.3,195,"a)",pos=4, cex=1.4)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)
text(0.3,195,"b)",pos=4, cex=1.4)
dev.off()
quartz_off_screen
3
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_figure1_v2.pdf", width = 10, height = 5)
par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^2*')'),ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)
text(0.3,690,"a)",pos=4, cex=1.4)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)
text(0.3,690,"b)",pos=4, cex=1.4)
dev.off()
quartz_off_screen
3
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/all_v3.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
legend(x="topleft", legend=c("human","lightning"),col="black",lty = c(1,2),pch=1,bty = "n",lwd = 2,cex=1.5, pt.cex = 1)
dev.off()
quartz_off_screen
3
t.test(log(pro1$growth_km),log(pro1h$growth_km))
Welch Two Sample t-test
data: log(pro1$growth_km) and log(pro1h$growth_km)
t = 6.5083, df = 169.36, p-value = 8.265e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
1.492386 2.791866
sample estimates:
mean of x mean of y
1.3226886 -0.8194377
t.test(log(pro2$growth_km),log(pro2h$growth_km))
Welch Two Sample t-test
data: log(pro2$growth_km) and log(pro2h$growth_km)
t = 6.9843, df = 134.4, p-value = 1.199e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
1.369607 2.451706
sample estimates:
mean of x mean of y
2.7559564 0.8453003
t.test(log(pro3$growth_km),log(pro3h$growth_km))
Welch Two Sample t-test
data: log(pro3$growth_km) and log(pro3h$growth_km)
t = 5.3694, df = 137.17, p-value = 3.286e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.9293526 2.0129234
sample estimates:
mean of x mean of y
3.186372 1.715234
t.test(log(pro4$growth_km),log(pro4h$growth_km))
Welch Two Sample t-test
data: log(pro4$growth_km) and log(pro4h$growth_km)
t = 4.6217, df = 122.92, p-value = 9.467e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.6909157 1.7261008
sample estimates:
mean of x mean of y
3.513204 2.304696
t.test(log(pro5$growth_km),log(pro5h$growth_km))
Welch Two Sample t-test
data: log(pro5$growth_km) and log(pro5h$growth_km)
t = 4.9328, df = 103.98, p-value = 3.089e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.9170094 2.1499736
sample estimates:
mean of x mean of y
3.895266 2.361774
##################3 for western cordillera ecoregion ##################
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2,cex.axis=1.4)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
null device
1
for mediteranean california
pro1 = res[res$fire_day == 1 & res$human == 1 & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & (res$eco1 == 11),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & (res$eco1 == 11),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="",cex.axis=1.4, xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
for difference in autumn
pro1 = res[res$fire_day == 1 & res$human == 1 & res$month >8,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month >8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month >8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month >8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month >8,]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & res$month >8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month >8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month >8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month >8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month >8,]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/autumn_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,250),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
for difference in summer
pro1 = res[res$fire_day == 1 & res$human == 1 & res$month <=8 ,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month <=8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month <=8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month <=8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month <=8,]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & res$month <=8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month <=8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month <=8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month <=8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month <=8,]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/summer_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
for difference in summer in western cordillera
pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
t.test(log(pro1$growth_km),log(pro1h$growth_km))
Welch Two Sample t-test
data: log(pro1$growth_km) and log(pro1h$growth_km)
t = 1.038, df = 53.703, p-value = 0.3039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4317907 1.3585976
sample estimates:
mean of x mean of y
-0.2855046 -0.7489080
t.test(log(pro2$growth_km),log(pro2h$growth_km))
Welch Two Sample t-test
data: log(pro2$growth_km) and log(pro2h$growth_km)
t = 3.6912, df = 82.447, p-value = 0.0003995
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.555548 1.854073
sample estimates:
mean of x mean of y
2.058908 0.854097
t.test(log(pro3$growth_km),log(pro3h$growth_km))
Welch Two Sample t-test
data: log(pro3$growth_km) and log(pro3h$growth_km)
t = 2.2426, df = 64.873, p-value = 0.02835
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.07573427 1.30867147
sample estimates:
mean of x mean of y
2.535161 1.842958
t.test(log(pro4$growth_km),log(pro4h$growth_km))
Welch Two Sample t-test
data: log(pro4$growth_km) and log(pro4h$growth_km)
t = 2.3845, df = 62.683, p-value = 0.02014
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1141261 1.2961201
sample estimates:
mean of x mean of y
2.976315 2.271192
t.test(log(pro5$growth_km),log(pro5h$growth_km))
Welch Two Sample t-test
data: log(pro5$growth_km) and log(pro5h$growth_km)
t = 2.3428, df = 40.216, p-value = 0.02417
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1060498 1.4368888
sample estimates:
mean of x mean of y
3.327351 2.555881
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_summer_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="",cex.axis=1.4, xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
for difference in autumn in western cordillera
pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_autumn_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,250),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
for difference in summer in meditereanean
pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
#t.test(log(pro1$growth_km),log(pro1h$growth_km))
#t.test(log(pro2$growth_km),log(pro2h$growth_km))
#t.test(log(pro3$growth_km),log(pro3h$growth_km))
#t.test(log(pro4$growth_km),log(pro4h$growth_km))
#t.test(log(pro5$growth_km),log(pro5h$growth_km))
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_summer_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
` ################## for difference in autumn in mediteranean ##################
pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)
me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)
points(me1,type="o",lty = 2, lwd = 2)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_autumn_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))
plot(me,type="o", ylim=c(0,400),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1,cex.axis=1.4, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off()
quartz_off_screen
3
how many days does it take to reach 75% burnt area
res75 = res[res$per_ba > 0.75,]
#peak_day = as.data.frame(aggregate(res75$fire_day, by = list(res75$firename,res75$cause), min))
#peak_day=subset(peak_day,x < 55)
#hi = hist(peak_day$x,prob =F, breaks= c(0:54), xlim=c(0,55), ylab="number of fires", xlab="days", cex.lab=1.4,cex.axis=1.3)
out1 = subset(res75,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 )
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak=as.data.frame(aggregate(res75$fire_day, by = list(res75$firename), min))
quantile(peak_day1$x,0.50,type=3)
50%
10
quantile(peak_day2$x,0.50,type=3)
50%
3
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/time_to_reach75.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
######## mediteranean ###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11)
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,15))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
######## north cal ###########
out1 = subset(res75,cause == 1 & (res75$eco1 == 6 | res75$eco1 == 7)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & (res75$eco1 == 6 | res75$eco1 == 7))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,15))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
######## mediteranean SUMMER ###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11 & ( res75$month <=8 & res75$month >5 )) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11 & ( res75$month <=8 & res75$month >5 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med_summer.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,10))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
######## north cal summer ###########
out1 = subset(res75,cause == 1 & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month <=8 & res75$month >5 )) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month <=8 & res75$month >5 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north_summer.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,12))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
######## mediteranean autumn###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11 & ( res75$month >8)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11 & ( res75$month >8 ))
#peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
#peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
#hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(0,hist.b$counts)
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med_autumn.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,5))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
######## north cal autumn ###########
out1 = subset(res75,cause == 1 & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month >8)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month >8 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north_autumn.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,5))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
2
res=res[res$ros_km>0,]
res_f = res[res$max_land == 1,]
res_p = res[res$max_land != 1,]
summary(lm(log(res$mean_frp)~log(res$ros_km)))
Call:
lm(formula = log(res$mean_frp) ~ log(res$ros_km))
Residuals:
Min 1Q Median 3Q Max
-2.5266 -0.5420 -0.0497 0.4990 4.5702
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.74609 0.01869 93.41 <2e-16 ***
log(res$ros_km) 0.35186 0.01321 26.63 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8362 on 2249 degrees of freedom
(547 observations deleted due to missingness)
Multiple R-squared: 0.2397, Adjusted R-squared: 0.2394
F-statistic: 709.1 on 1 and 2249 DF, p-value: < 2.2e-16
#just show the plot here
plot(res_f$mean_frp~res_f$ros_km,log="xy",xlim=c(0.005,30),ylim=c(0.1,180),xaxt="n",ylab="mean FRP (MW)",xlab="Rate-of-Spread (km/day)", cex.lab=1.4,cex.axis = 1.3,col="darkgreen")
marks=c(0.01,0.1,1,10)
marks1=c(0.1,0.5,5,50)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/fig_FRP_ros_v3.pdf", width = 5, height = 5)
plot(res_f$mean_frp~res_f$ros_km,log="xy",xlim=c(0.005,30),ylim=c(0.1,180),xaxt="n",yaxt="n",ylab="mean FRP (MW)",xlab=expression('Rate-of-Spread (km d'^-1*')'), cex.lab=1.4,cex.axis = 1.3,col="darkgreen")
points(res_p$mean_frp~res_p$ros_km,col="orange")
axis(1,at=marks,labels=marks,cex.axis=1.4 )
axis(2,at=marks1,labels=marks1,cex.axis=1.4 )
legend( x="topleft",legend=c("Forest","Grass & shrub"),col=c("darkgreen","orange"),cex=1.2,pch=1,bty = "n")
dev.off()
quartz_off_screen
3
FRP vs tree mortality
summary(lm(log(res$mean_frp)~res$mean_bas))
Call:
lm(formula = log(res$mean_frp) ~ res$mean_bas)
Residuals:
Min 1Q Median 3Q Max
-3.3250 -0.5824 0.0428 0.5755 2.7943
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.853006 0.039957 21.35 <2e-16 ***
res$mean_bas 0.019447 0.001025 18.97 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8574 on 1727 degrees of freedom
(1069 observations deleted due to missingness)
Multiple R-squared: 0.1725, Adjusted R-squared: 0.172
F-statistic: 359.9 on 1 and 1727 DF, p-value: < 2.2e-16
summary(lm(data_forest$mean_BA_red~(log10(data_forest$ros_km)+I((log10(data_forest$ros_km))^2))))
Call:
lm(formula = data_forest$mean_BA_red ~ (log10(data_forest$ros_km) +
I((log10(data_forest$ros_km))^2)))
Residuals:
Min 1Q Median 3Q Max
-40.210 -11.114 -3.598 8.209 82.558
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 27.0251 0.5005 54.00 <2e-16 ***
log10(data_forest$ros_km) 25.2022 1.1296 22.31 <2e-16 ***
I((log10(data_forest$ros_km))^2) 8.3470 0.8597 9.71 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.51 on 1426 degrees of freedom
Multiple R-squared: 0.2908, Adjusted R-squared: 0.2898
F-statistic: 292.4 on 2 and 1426 DF, p-value: < 2.2e-16
difference in fire size for first 5 days across california and both ecosystems
res$ros1 = res$max_ros+1
out1 = subset(res,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )
hum1 = out2[out2$fire_day ==1,]
hum2 = out2[out2$fire_day ==2,]
hum3 = out2[out2$fire_day ==3,]
hum4 = out2[out2$fire_day ==4,]
hum5 = out2[out2$fire_day ==5,]
lig1 = out1[out1$fire_day ==1,]
lig2 = out1[out1$fire_day ==2,]
lig3 = out1[out1$fire_day ==3,]
lig4 = out1[out1$fire_day ==4,]
lig5 = out1[out1$fire_day ==5,]
mean(hum1$growth)
mean(hum2$growth)
mean(hum3$growth)
mean(hum4$growth)
mean(hum5$growth)
mean(lig1$growth)
mean(lig2$growth)
mean(lig3$growth)
mean(lig4$growth)
mean(lig5$growth)
t.test(log10(hum1$growth),log10(lig1$growth))
t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
t.test(log10(hum4$growth),log10(lig4$growth))
t.test(log10(hum5$growth),log10(lig5$growth))
hum1 = out2[out2$fire_day ==1 & out2$eco1==6,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==6,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==6,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==6,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==6,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==6,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==6,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==6,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==6,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==6,]
mean(hum1$growth)
mean(hum2$growth)
mean(hum3$growth)
mean(hum4$growth)
mean(hum5$growth)
mean(lig1$growth, na.rm=T)
mean(lig2$growth, na.rm=T)
mean(lig3$growth, na.rm=T)
mean(lig4$growth, na.rm=T)
mean(lig5$growth, na.rm=T)
t.test(log10(hum1$growth),log10(lig1$growth))
t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
t.test(log10(hum4$growth),log10(lig4$growth))
t.test(log10(hum5$growth),log10(lig5$growth))
hum1 = out2[out2$fire_day ==1 & out2$eco1==11,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==11,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==11,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==11,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==11,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==11,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==11,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==11,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==11,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==11,]
mean(hum1$growth)
mean(hum2$growth)
mean(hum3$growth)
mean(hum4$growth)
mean(hum5$growth)
mean(lig1$growth, na.rm=T)
mean(lig2$growth, na.rm=T)
mean(lig3$growth, na.rm=T)
mean(lig4$growth, na.rm=T)
mean(lig5$growth, na.rm=T)
#t.test(log10(hum1$growth),log10(lig1$growth))
#t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
t.test(log10(hum4$growth),log10(lig4$growth))
t.test(log10(hum5$growth),log10(lig5$growth))
out1 = subset(res,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )
mean(out1$ros_km,na.rm=T)
[1] 0.8405408
mean(out2$ros_km,na.rm=T)
[1] 1.82709
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
out1 = subset(res,cause == 1 & eco1 == 11) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11)
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
############ autumn northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 9)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 9))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
############ autumn mediterean california
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 9)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 9))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
############ summer northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
############ summer mediterean california
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 5 & month<10)) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 5 & month<10))
hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
out1 = subset(res,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )
mean(out1$ros_km,na.rm=T)
[1] 0.8405408
mean(out2$ros_km,na.rm=T)
[1] 1.82709
#0,0.25,0.5,1,2,3,5,7,10,20,30
hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure2_freq_v1.pdf", width = 6, height = 5)
fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="Fire days (%)",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
3
###burnt area by rate-of-spread
data_test = data_s1[data_s1$max_land == 1,]
data_test1 = data_s1[data_s1$L1CODE == 6 |data_s1$L1CODE == 7 ,]
data_test2 = data_s1[data_s1$L1CODE == 11,]
data_test1 = data_s1[data_s1$human == 1 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]
data_test2 = data_s1[data_s1$human == 2 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]
data_test1 = data_s1[data_s1$human == 1 ,]
data_test2 = data_s1[data_s1$human == 2 ,]
fast = data_test1[data_test1$ros_km > 3.64,]
sum(fast$size,na.rm=T)/sum(data_test1$size,na.rm=T)
[1] 0.4992649
sum((data_test1$size*data_test1$ros_km),na.rm=T)/sum(data_test1$size,na.rm=T)
[1] 5.610449
fast = data_test2[data_test2$ros_km > 2.2,]
sum(fast$size,na.rm=T)/sum(data_test2$size,na.rm=T)
[1] 0.499989
sum((data_test2$size*data_test2$ros_km),na.rm=T)/sum(data_test2$size,na.rm=T)
[1] 3.532566
fast = data_s1[data_s1$ros_km > 1,]
fast_hum = fast[fast$human == 1,]
print("% BA and % of fire days fast fires > 1km/day")
[1] "% BA and % of fire days fast fires > 1km/day"
sum(fast$size)/sum(data_s1$size)
[1] 0.8359887
length(fast$size)/length(data_s1$size)
[1] 0.3371951
print("% BA fast fires due to human ignition % of fire days human caused fast fires > 1km/day")
[1] "% BA fast fires due to human ignition % of fire days human caused fast fires > 1km/day"
sum(fast_hum$size, na.rm=T)/sum(fast$size)
[1] 0.4634998
length(fast_hum$size)/length(fast$size)
[1] 0.4195298
all_min1 = data_s1[data_s1$nr_day != 1,] # remove first fire spread day from statistics
quan = quantile(data_s1$ros_km,0.9)
fast = data_s1[data_s1$ros_km > quan,]
slow = data_s1[data_s1$ros_km < quan,]
fast_hum = fast[fast$human == 1,]
print("fastest 10% fires cause xxx% of BA")
[1] "fastest 10% fires cause xxx% of BA"
sum(fast$size)/sum(all_min1$size)
[1] 0.5493523
length(fast$size)/length(all_min1$size)
[1] 0.1
print("mean tree mortality weighted by BA and just mean")
[1] "mean tree mortality weighted by BA and just mean"
sum((data_s1$mean_BA_red*data_s1$size))/(sum(data_s1$size))
[1] 49.03327
mean(data_s1$mean_BA_red)
[1] 25.64874
print("% BA due to human fires amoung fastest 10% fire days")
[1] "% BA due to human fires amoung fastest 10% fire days"
sum(fast_hum$size, na.rm=T)/sum(fast$size)
[1] 0.5143879
print("% fire number due to human fires amoung fastest 10% fire days")
[1] "% fire number due to human fires amoung fastest 10% fire days"
length(fast_hum$size)/length(fast$size)
[1] 0.5060976
print("% tree mortality fast fires weighthed and not")
[1] "% tree mortality fast fires weighthed and not"
sum((fast$mean_BA_red*fast$size))/(sum(fast$size))
[1] 60.02651
mean(fast$mean_BA_red)
[1] 52.79125
print("% tree mortality slow fires weighthed and not")
[1] "% tree mortality slow fires weighthed and not"
sum((slow$mean_BA_red*slow$size))/(sum(slow$size))
[1] 35.63221
mean(slow$mean_BA_red)
[1] 22.63291
print("tree mortality <0.5km and >2km")
[1] "tree mortality <0.5km and >2km"
fast1 = data_s1[data_s1$ros_km > 2,]
slow1 = data_s1[data_s1$ros_km < 0.5,]
mean(fast1$mean_BA_red,omit.na=T )
[1] 48.30526
mean(slow1$mean_BA_red,omit.na=T )
[1] 15.29343
# plot BA per rate of spread
out1 = subset(data_s1,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(data_s1,cause !=1 & cause != 14 )
breaks =c(0,0.5,1,2,3,5,7,10,20,30)
tt=0
pp=0
for (i in 1:9){
kr = out1[out1$ros_km >= breaks[i] & out1$ros_km < breaks[i+1],]
kp = out2[out2$ros_km >= breaks[i] & out2$ros_km < breaks[i+1],]
tt[i]=sum(kr$size, na.rm=T)/1000000
pp[i]= sum(kp$size, na.rm=T)/1000000
}
sum(tt)
[1] 6781.34
sum(pp)
[1] 5612.935
fg = rbind(tt,pp)
par(mar=c(4, 5, 2,0.1))
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab=expression('burnt area (km'^2*')'),ylim=c(0,2000),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels= breaks ,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/BA_per_ros_hist.tif", width = 6, height = 5, units = 'in', res = 300)
par(mar=c(4, 5, 2,0.1))
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab=expression('burnt area (km'^2*')'),ylim=c(0,2000),cex.lab=1.4,cex.axis = 1.3)
#axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels= breaks ,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
3
fg1=rbind((tt/sum(tt))*100,(pp/sum(pp))*100)
pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/BA_per_ros_hist_percentage.pdf", width = 6, height = 5)
#par(mar=c(1, 1, 1,1))
fr = barplot(fg1, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab= "Burned area (%)",ylim=c(0,30),cex.lab=1.4,cex.axis = 1.3)
#axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels= breaks ,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()
quartz_off_screen
3
20 fastest fires
res1 = res[res$ros_km > 10 & !is.na(res$ros_km),]
Error in res$ros_km : object of type 'closure' is not subsettable
are ROS the same for light & human under the same conditions
out1 = subset(res,cause == 1 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )
plot(out1$vpd,log(out1$ros_km))
points(out2$vpd,log(out2$ros_km),col="red")
summary(lm(out1$vpd~log(out1$ros_km+1),na.omit=T))
summary(lm(out2$vpd~log(out2$ros_km+1)))
analysis of the first day
load data
daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/all_ignitions_V3.txt",header=T)
res=as.data.frame(daily_res)
res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$total_area =as.numeric(as.character(res$total_area))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))
res$biomass =as.numeric(as.character(res$biomass))
res = res[-1,]
res$human[res$cause ==1] =1
res$human[res$cause !=1 & res$cause !=14] =0
analysis
out1 = res[res$cause !=1 & res$cause != 14,]
out2 = res[res$cause ==1,]
length(out1$bi)
[1] 1196
length(out2$bi)
[1] 486
out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1) #1=lightning; 14=unknown; 7=arson
length(out1$bi)
[1] 335
length(out2$bi)
[1] 375
out1 = subset(res,eco1 == 11& res$cause !=1 & res$cause != 14)
out2 = subset(res,eco1 == 11 & res$cause ==1)
length(out1$bi)
[1] 784
length(out2$bi)
[1] 45
out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14 & res$mont > 5 & res$mont < 10 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1& res$mont > 5 & res$mont < 10) #1=lightning; 14=unknown; 7=arson
length(out1$bi)
[1] 248
length(out2$bi)
[1] 352
out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14 & res$mont < 6 & res$mont > 9 ) #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1& res$mont < 6 & res$mont > 9) #1=lightning; 14=unknown; 7=arson
length(out1$bi)
[1] 0
length(out2$bi)
[1] 0
t.test(out1$bi,out2$bi)
Error in t.test.default(out1$bi, out2$bi) : not enough 'x' observations
analysis of MTBS temporal trend
library(raster)
library(rgeos)
mtbs_dir="/Users/stijnhantson/Documents/data/MTBS/DATA/"
med=shapefile("/Users/stijnhantson/Documents/data/veg_california/med_cal.shp")
north=shapefile("/Users/stijnhantson/Documents/data/veg_california/north_cal.shp")
year=1984
k=0
dnbr_all_med=0
rdnbr_all_med=0
dnbr_all_nor=0
rdnbr_all_nor=0
len_dnbr_med=0
len_rdnbr_med=0
len_dnbr_nor=0
len_rdnbr_nor=0
for (year in 1984:2017){
print(year)
k=k+1
mtbs_dir_year = paste(mtbs_dir,year,"/",sep="")
bndy_list = list.files(mtbs_dir_year, pattern = "burn_bndy.shp$", recursive = T, full.names=T)
shapefile_list <- lapply(bndy_list, shapefile)
fires <- do.call(rbind, shapefile_list)
fires=gUnaryUnion(fires)
fire_north = intersect(fires,north)
fire_med = intersect(fires,med)
dnbr = raster(paste(mtbs_dir,year,"_dnbr.tif",sep=""))
dnbr[dnbr < -2000] <- NA
rdnbr = dnbr = raster(paste(mtbs_dir,year,"_rdnbr.tif",sep=""))
rdnbr[rdnbr < -2000] <- NA
dnbr_ext_med = extract(dnbr,fire_med)
dnbr_ext_nor = extract(dnbr,fire_north)
rdnbr_ext_med = extract(rdnbr,fire_med)
rdnbr_ext_nor = extract(rdnbr,fire_north)
dnbr_all_med[k]=mean(unlist(dnbr_ext_med),na.rm=T)
rdnbr_all_med[k]=mean(unlist(rdnbr_ext_med),na.rm=T)
dnbr_all_nor[k]=mean(unlist(dnbr_ext_nor),na.rm=T)
rdnbr_all_nor[k]=mean(unlist(rdnbr_ext_nor),na.rm=T)
len_dnbr_med[k] = length(unlist(dnbr_ext_med))
len_rdnbr_med[k] = length(unlist(rdnbr_ext_med))
len_dnbr_nor[k] = length(unlist(dnbr_ext_nor))
len_rdnbr_nor[k] = length(unlist(rdnbr_ext_nor))
removeTmpFiles(0)
gc()
}
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
```